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This paper develops goodness of fit statistics that can be used to 
formally assess Markov random field models for spatial data, when 
the model distributions are discrete or continuous and potentially 
parametric. Test statistics are formed from generalized spatial resid- 
uals which are collected over groups of nonneighboring spatial ob- 
servations, called concliques. Under a hypothesized Markov model 
structure, spatial residuals within each conclique are shown to be in- 
dependent and identically distributed as uniform variables. The infor- 
mation from a series of concliques can be then pooled into goodness 
of fit statistics. Under some conditions, large sample distributions 
of these statistics are explicitly derived for testing both simple and 
composite hypotheses, where the latter involves additional paramet- 
ric estimation steps. The distributional results are verified through 
simulation, and a data example illustrates the method for model as- 
sessment. 



1. Introduction. Conditionally specified models formulated on the basis 
of an underlying Markov random field (MRF) are an attractive alternative 
to continuous random field specification for the analysis of problems that in- 
volve spatial dependence structures. By far the most common of such models 
are those formulated using a conditional Gaussian distribution (e.g., [42]), 
but models may also be constructed using a number of other conditional 
distributions such as a beta [23, 31], binary [8], Poisson [4] or Winsorized 
Poisson [29], and general specifications are available for many exponential 
families [2, 31]. 
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In an applied spatial setting, we assume that observations are available 
at a finite set of geo-referenced locations {sj:i = 1, ...,N}, and to these 
locations we assign the random variables {Y(s,) : i = 1, . . . , N}. In general, 
locations are arbitrarily indexed in d-dimensional real space. A MRF is 
typically constructed by specifying for each location Sj a neighborhood, con- 
sisting of other locations on which the full conditional distribution of Y(sj) 
will be functionally dependent. Let the conditional cumulative distribution 
function (c.d.f.) of Y(sj) given {Y(sj) = y(sj) -j^i} be denoted as Fi and 
define Mi = {sj ^ Sj, and Fi depends functionally on y(sj)}. Also define 

y(M) = {y( 

Sj) :sj eMi}. The Markov assumption implies that 

(1.1) F i (-\{y( Sj ): Sj ^s i })=F i ^{y(s j ): Sj eM i })=F i (-\y(M i )). 

A model is formulated by specifying, for each i = 1, . . . , N, a conditional 
c.d.f. in (1.1). Conditions necessary for a set of such conditionals to corre- 
spond to a joint distribution for {Y(si), . . . , Y(sjv)} are given by Arnold, 
Castillo and Sarabia [2] and a constructive process with useful conditions 
sufficient for existence of a joint are laid out in Kaiser and Cressie [30]. Mod- 
els may be constructed for both discrete and continuous random variables, 
on regular or irregular lattices, with or without an equal number of neigh- 
bors for each location (including Mi = for some locations) and possibly 
including information from spatial covariates. The construction of models 
for applications is thus very flexible. 

A number of our results and, in particular, Theorem 2.1 to follow, can 
be generalized to some of the variable situations just described, but it will 
be beneficial for developing theoretical results to define a setting that is 
broad but highly structured. We desire a spatial process defined on grid 
nodes of the d-dimensional integer lattice Z rf , where Z = {0, ±1, ±2, . . .}. We 
stipulate a number of restrictions for this process that, while not capable of 
covering all of the finite-dimensional models mentioned previously, is flexible 
enough to be meaningful in many applied situations. We formally consider 
specifying an MRF model for a spatial process Y = {Y(s) : s € Z d }, rather 
than a model (1.1) developed with respect to a finite collection of (possi- 
bly nonlattice) data sites {Y(sj) : i = 1, . . . , N}. To this end, assume that 
for any s € Z rf neighborhoods can be constructed using a standard template 
M C Z d \{0} as M(s) = s + M, with \M\ < oo denoting the size of M. Some 
examples of A4 are given in the next section. We then assume that the pro- 
cess Y has a stationary distribution function F(-\-) such that, for any s € Z d , 
the conditional c.d.f. of Y(s) given all remaining variables {Y(t):t € 
t ^ s} can be written as 

(1.2) F(-\{Y(t) :t G Z d ,t + s}) = F(-\{Y(t) :t € M(s)}) 
under a Markov assumption. 
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Given a hypothesized or estimated model, our concern is how one might 
conduct a goodness of fit (GOF) procedure, either through informal diag- 
nostics or by using formal probability results that lead to a GOF test. The 
approach we propose here may be viewed within either the context of a pure 
GOF test to address the question of whether a (possibly fitted) model pro- 
vides an adequate description of observed data. This is an issue of model 
assessment and different from model selection, which has been considered, 
for example, with penalized pseudo-likelihood for parametric MRF mod- 
els; cf. [11, 21, 26]. Additionally, while other GOF tests may be possible 
for certain joint model specifications (e.g., a frequency-domain approach for 
Gaussian processes; cf. [1]), we focus solely on conditional model specifica- 
tions. The GOF variates introduced in the next section may be used as either 
diagnostic quantities or as the basis for a formal GOF test as presented in 
Section 3. 

The remainder of this article is organized as follows. In Section 2 we in- 
troduce the concept of a conclique and derive GOF variates that form the 
basis of our approach, using an adaptation of a multivariate probability 
integral transform (PIT). Section 3 develops a formal methodology for com- 
bining these variates over concliques to create GOF tests of Markov models 
under both simple and composite hypotheses. These tests are omnibus in 
the sense that they assess the hypothesized model in total, including the 
neighborhood structure selected, specification of dependence as isotropic or 
directional, and the form of the modeled conditional distributions. Theoret- 
ical results are presented in Section 4 that establish the limiting sampling 
distributions of GOF tests under the null hypothesis. Section 5 describes 
a numerical study to support the theoretical findings. Section 6 provides an 
application of the GOF tests in model assessment for agricultural trials. Sec- 
tion 7 contains concluding remarks and discussions on extensions. Section 8 
provides a proof of the foundational conclique result (Theorem 2.1), and all 
other proofs regarding the asymptotic distribution of GOF test statistics 
appear in supplementary material [32]. 

2. Generalized spatial residuals. In this section we derive the basic quan- 
tities that form the basis for our GOF procedures. We consider these quan- 
tities to be a type of generalized residuals because they fit within the frame- 
work suggested by Cox and Snell [9]. In particular, these generalized spatial 
residuals will be derived using an extended version of Rosenblatt's [41] mul- 
tivariate PIT combined with a partitioning of spatial locations into sets such 
that the residuals within each set constitute a random sample from a uni- 
form distribution on the unit interval, under the true model. As discussed 
by Brockwell [7] and Czado et al. [12], the PIT formulation allows arbitrary 
model distributions to be considered in assessing GOF, rather than simply 
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continuous ones. Similar transformations, with subsequent formal or infor- 
mal checks for uniformity, have been important in evaluating the GOF of, 
and the quality of predictive forecasts from, various models for time series; 
cf. [13, 15, 16, 19, 24, 27]. 

2.1. Concliques. Before providing the transform that defines our gener- 
alized spatial residuals, it is necessary to develop a method for partitioning 
the total set of spatial locations at which observations are available into sub- 
sets with certain properties. We call such sets concliques because they are 
defined as the converse of what are called cliques by Hammersley and Clif- 
ford [22] . In the case of regular lattices with neighborhoods defined using ei- 
ther four-nearest or eight-nearest neighbor structures, concliques correspond 
exactly to the so-called coding sets of Besag [4], which were suggested for 
use in forming conditional likelihoods for estimation. The key property of 
concliques, however, allows construction of such sets in more general settings 
including irregular lattices and hence the new name. 

As defined in [22], a clique is a set of locations such that each location in 
the set is a neighbor of every other location in the set. Similar terminology 
exists in graph theory, where a subset of graph vertices (e.g., locations) form 
a clique if every two vertices in the subset are connected by an edge [45]. 
We define a conclique as a set of locations such that no location in the set is 
a neighbor of any other location in the set. Any two members of a conclique 
may share common neighbors, but they cannot be neighbors themselves. 
Additionally, every set of a single location can be treated as both a clique or 
conclique. In the parlance of graphs, the analog of a conclique is a so-called 
"independent set," defined by a set of vertices in which no two vertices share 
an edge. This particular graph terminology conflicts with the probabilistic 
notion of independence, while a "conclique" truly represents a conditionally 
independent set of locations in a MRF model. 

While the result of the next subsection holds for any collection of con- 
cliques, in practice what is desired is a collection of concliques that suit- 
ably partition all observed locations. To achieve this under the process 
model (1.2), we identify a collection of concliques {Cj :j = l,...,q} that par- 
tition the entire grid Z rf . We define a collection of concliques to be a minimal 
conclique cover if it contains the smallest number of concliques needed to 
partition the set of all locations. In graph theory, this concept is related to 
determining the smallest (or chromatic) number of colors needed to color 
a graph (with no two edge-connected vertices sharing the same color) or, 
equivalently, the smallest number of independent sets needed to partition 
graph vertices [25]. In practice, identifying a minimal conclique cover is 
valuable since our procedure produces one test statistic for each conclique 
in a collection, and those statistics must then be combined into one overall 
value for a formal GOF test. 
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Example 2.1 (A 4-nearest neighbor model on Z 2 ). Here, let s = (it, v)' € 
Z 2 for a horizontal coordinate u and a vertical coordinate v. The neighbor- 
hood structure of a 4-nearest neighbor model is produced with the template 
M = {(-1,0)', (1,0)', (0,1)', (0,-1)'}, so that M(s) for a given location s 
and neighbors * is as shown in the following figure: 

* 

* s * 

* 

In this case, the minimal conclique cover contains two members, C\ and C2, 
with elements denoted by l's and 2's, respectively, as shown below. 
Minimal conclique collection for a 4-nearest neighbor model: 

12 12 12 12 12 
2 12 12 12 12 1 
12 12 12 12 12 
2 12 12 12 12 1 
12 12 12 12 12 

Example 2.2 (An 8-nearest neighbor model on Z 2 ). As in the previ- 
ous example, let s = (it, f)' but take M. = {(u,v)' :max{|it|, \ v\} = 1}. The 
neighborhood structure of an 8-nearest neighbor model is then shown in the 
following figure for a location s € Z 2 and neighbors *: 

# # # 

* s * 

# # # 

For the 8-nearest neighbor model, there are four concliques in the mini- 
mal cover, C\, . . . ,£4, with elements denoted by l's, 2's, 3's and 4's in the 
following figure, respectively. 

Minimal conclique cover for an 8-nearest neighbor model: 

12 12 12 12 12 
3434343434 
12 12 12 12 12 
3434343434 
12 12 12 12 12 

2.2. Defining generalized spatial residuals. Let {^4(s) : s € Z d } denote a col- 
lection of independent and identically distributed (i.i.d.) random variables, 
which are Uniform (0, 1) and also independent of the spatial process Y. For 
any s 6 Z rf , we then define a random generalized spatial residual as 

U(b) = (1 - A(s)) ■ F(Y(s)\{Y (t) : t G M(s)}) 
+ A(s)-F-(Y(s)\{Y(t):teM(s)}), 
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where F(-\-) denotes the (stationary) c.d.f. from (1.2), and F (-|-) denotes 
the left limit of the c.d.f., that is, F~(y\{Y(t) :t G W(s)}) = P(Y(s) < y\ 
{Y(t):t GJV(s)}), y € M. This residual applies the notion of a random- 
ized PIT [7], allowing for a noncontinuous c.d.f. .F(-|-) to be considered. 
When F(-\-) is continuous, the spatial residual reduces to a PIT U(s) = 
F(Y(s)\{Y(t):t e7V(s)}) in Rosenblatt's [41] format. Given that a collec- 
tion of concliques is available for a particular situation, the fundamental 
result that serves as the basis for our GOF procedures is as follows. 

Theorem 2.1. Let the spatial process {Y(s):s€Z a! } have conditional 
distribution functions as in (1.2), and let {Cj : j = 1, . . . ,q} be a collection of 
concliques that partition the integer grid 7h d . Then for any j = 1, .. . ,q, the 
variables {C/(s):s€Cj} given by (2.1) are i.i.d. Uniform (0,1) variables. 

Typically, the conditional c.d.f. F(-\-) of expression (1.2) will be a pa- 
rameterized function, and we now write this as Fg(-\-) to emphasize the 
parametrization. Let 0q denote the true value of the parameter. In an appli- 
cation we have available a set of observations taken to represent realizations 
of the random variables {Y(sj) :i = 1, .. . ,N}. Theorem 2.1 indicates that if 
we compute generalized spatial residuals as, in the notation of (2.1), 



then within any conclique Cj these variables should behave as a random 
sample from a uniform distribution on the unit interval. If we use a minimal 
conclique cover having q members, then we will have q sets of residuals, each 
of which should behave as a random sample from a uniform distribution. 
These sets of residuals will not, however, be independent, so we will not 
have a total collection that behaves as q independent random samples. 

In practice we will usually also replace the parameter 6 with an estimate 9 
computed on the basis of the observations so that, technically, the values 
within any conclique will not actually be independent either. We expect, 
however, that if the model is appropriate, then these residuals will exhibit 
approximately the same behavior as independent uniform variates, in the 
same way that ordinary residuals from a linear regression model with normal 
errors behave as an approximate random sample of normal variates, despite 
the fact that they cannot technically represent such a sample. 

A basic diagnostic plot can be constructed by plotting the empirical dis- 
tribution function of each set of residuals {tt(s,) : Sj £ Cj}, j = 1, . . . , q, and 
examining them for departures from a standard uniform distribution func- 
tion. See, for instance, Gneiting et al. [19], Section 3.1, for a summary of 



(2.2) 



U( Si ) = (1 - Afa)) ■ FQ (y(si)\{y(t) :t eJV(sj)}) 
+ A( Si )-F 9 -(y( Si )\{y(t):teM( Si )}), 



Si e Cj 
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graphical approaches for exploring uniformity in PIT values. Tests for unifor- 
mity may be used for individual sets of residuals to guide the decision about 
whether a given fitted model is adequate or to choose between two compet- 
ing (even nonnested) models. Such procedures do not constitute a formal 
GOF test, however, because there is no guarantee that results will agree 
across differing sets of residuals in a conclique cover. Formal procedures for 
combining evidence from the residual sets into one overall GOF test are 
presented in the next section. 

3. Methodology: Goodness of fit tests. 

3.1. General setting. Suppose that for a set of locations on the d-dimen- 
sional integer lattice {si, . . . , sjy} C we want to assess the GOF of a con- 
ditional model specification, based on a set of observed values {V(sj) :i = 
1,...,N}. We assume that the observed values are a partial realization of 
a class of process models defined on Z, d for which the conditional c.d.f. 
of y(s) given {Y(t) :t ^s} belongs to a class of parameterized conditional 
distribution functions, 

(3.1) T e = {F e (-\{Y(t) : t e JV(s)}) :9e8}, 

where 0CF, 1 < p < oo, is a parameter space, Af(s) = s + M and, analo- 
gously to (1.2), M <Z r L d \ {0}. Two testing problems fit into this framework, 
where the null hypothesis is simple and where it is composite. 

In the next subsections, we describe GOF tests for simple and compos- 
ite hypotheses based on the observations {Y(sj):i = 1, ...,N}, which are 
assumed to have arisen in the following way. Suppose that R C M. d denotes 
a sampling region within which N observations are obtained at a set of sam- 
pling locations Sn = R H Z rf = {si, . . . , sjv}- Define the interior of the set of 
sampling locations as Sffi = {s£ Sn :N(s) C <Stv}. Locations in this set are 
those sampling locations for which all neighbors are also sampling locations, 
allowing generalized spatial residuals to be computed for all s 6 Sffi , even if 
the physical sampling region R is irregular. Finally, let Cin, . . . ,C q N denote 
the conclique partition of Sjy* determined by CjN = Cj fl Sffi, j = 1, • • • , <?• In 
practice we will desire a minimal conclique cover but this is not necessary 
in what follows. 

3.2. Testing a simple null hypothesis. First consider the case of the sim- 
ple (S) null hypothesis in which the testing problem is given by, for some 
specified 8qG@, 

Hq(S): The data {Y(sj) :i = 1, . . . , N} represent a partial sample of 

the process model class (3.1) with 6 = 6q; 
#1 (5): Not H (S). 
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To construct test statistics appropriate for these hypotheses, we consider 
the generalized spatial residuals under Hq(S), 

U(s) = (1 - A(s)) ■ F eo (Y(s)\{Y(t) : t G M(s)}) 

3 ' 2 + A(a) ■ F e -(Y(s)\{Y(t) :t G M(s)}), s G S%\ 

Now define, for j = 1, . . . ,q, the (generalized residual) empirical distribution 
function over the jth conclique by 



G jN {u) = -±- £ l(U(s)<u), 



\CjN\ 

u G [0, 1]. Here and in the following, 1(A) denotes the indicator function of 
a statement A, where 1(A) = 1 if A is true and 1(A) = otherwise. Note that 
under Hq(S), E{GjN(u)} = u, u G [0, 1], as a result of Theorem 2.1. Hence, 
to assess the GOF of the model over the jth. conclique Cj , we consider the 
scaled deviations of the empirical distribution function from the Uniform 
(0, 1) distribution, 

(3.3) W jN (u) = N l ' 2 (G jN (u)-u), uG[0,l]. 

A number of GOF test statistics for testing Hq(S) may be obtained by 
combining the Wj^'s in different ways: 

(3.4) Tin = max sup |WjAr(tt)|, 

j=i,-,g ne [ ,i] 



/ 1 

(3.5) T 2N = - V sup \W jN (u)\ 

V«^ L «6[0,1] 

(3.6) T3N = max I / \WjN(u)\ r du 

J=h-,q\Jo 

(3.7) T m = ^f^^\W jN (u)\ r d^j 



1/2 



l/r 



where r G [l,oo) in (3.6) and (3.7). Note that Tin and T27V are obtained 
by combining conclique- wise Kolmogorov-Smirnov test statistics, while T^n 
and T47V are obtained by combining conclique-wise Cramer-von Mises test 
statistics. While our statistics are based exclusively on paired differences 
(e.g., GjN(u) — u, u G [0, 1]), other test statistics may be formulated to assess 
agreement between the empirical GjN and Uniform(0, 1) distributions, such 
as GOF tests based on (^-divergences studied in [24]. In Section 4, we provide 
asymptotic distributions for the empirical processes (3.3), which may be 
an ingredient for determining limit distributions of statistics based on (p- 
divergences; cf. Theorem 3.1 [24]. 
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3.3. Testing a composite null hypothesis. The composite (C) null hy- 
pothesis can be stated as 

Hq{C): The data {Y(sj) : i = 1, . . . , N} represent a partial sample of 

some member of the process model class (3.1) for an unknown 9; 

H X {C): Not H Q {C). 

Let 9 denote an estimator of 9 based on {Y(sj) : i = 1, . . . , N}. Since 9 is 
unknown, instead of the ?7(s)'s of (3.2), we work with an estimated version 
of the generalized spatial residuals, 

U(s) = (1 - A(s)) ■ F § (Y(s)\{Y(t) : t G Af(s)}) 

+ A(s) ■ F7 (Y (s)\{Y (t) :t G AA(s)}), s G S]f , 

where, as before, A/"(s) = s + A4. Note that if is a reasonable estimator 
of 9 and if is a smooth function of 0, then the C/(s)'s of (3.8) are 

approximately distributed as Uniform (0, 1). This suggests that we can base 
tests of Hq(C) versus H\(C) on the processes 

(3.9) W jN {u) = N l / 2 (G jN (u)-u), u€[0,l], 
for j = 1, ... j g, where 

= 17^7 E E ^( s ) ^ n )' « e [0, !]■ 

The test statistics for testing Hq(C) versus H\(C) are now given by 

(3.10) Tiat, . . . ,T 4 at, 

where TjN is obtained by replacing Wjjy in expressions (3.4)-(3.7) with WjN- 
In the next section, we describe the limit distributions of the test statistics 
under the null hypothesis. 

4. Asymptotic distributional results. 

4.1. Basic concliques. To formulate large sample distributional results 
for the GOF statistics, we shall assume that the concliques C\,...,C q used 
for these statistics can be "built up" from unions of structurally more basic 
concliques, say C*,...,C**, q* > q. For any given template M C 7L d \ {0} 
defining neighborhoods as Af(s) = s + M., s G Z d , we suppose such concliques 
are constructed as follows. 

Let ej G Z d denote a vector with 1 in the ith component and elsewhere, 
and define rrii = max{|e^s| :s G M} as the maximal absolute value of ith 
component over integer vectors in the neighborhood template s G M, i = 
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1, . . . ,d. Define a collection of sublattices as 

d 

(4.1) C* = { aj +As:seZ d }, j = l,...,q* = l[(m i + l), 

i=l 

where A = diag(mi + 1, . . . , + 1) is a positive diagonal matrix and 

&j GX = {(ai, . . . ,ad)' € Z d :0 < aj < m,i,i = l,...,d}, 

where a, ^ a k if C* / C* k . 

Proposition 4.1 shows that these sets provide a collection of "basic" con- 
cliques (or coding sets) since locations within the same sublattice C* are 
separated by directional distances A that prohibit neighbors within C* . Ad- 
ditionally, the proposition gives a simple rule for merging basic concliques Cj 
to create larger concliques Cj. In the following, write ±A4 = A4 U — M, and 
define HsHoo = maxi<j<d |sj| for s = (s\, . . . , s^)' € Z d . 

Proposition 4.1. Under the process assumptions of Theorem 2.1 and 
for any neighborhood specified by a finite subset M. C Z rf \ {0}: 

(a) sets C\, . . . ,C*« of form (4-1) are concliques that partition Z d ; 

(b) if&i, . . . , aj, aj+i € I, i > 1, siic/i i/iai C = U}=i ^ s a conclique, then 
C U C 4 * +1 is a conclique if and only if 

&j — aj + i + As ^ ±Ai for all s S Z d , ||s||oo < 1, and any j = 1, . . . ,i. 

In addition to providing a systematic approach for building concliques, 
the purpose of this basic conclique representation is to allow the covariance 
structure of the limiting Gaussian process of the conclique-wise empirical 
processes [cf. (3.3)] to be written explicitly and to simplify the distributional 
results to follow (as basic concliques C* above have a uniform structure and 
are translates of one another). With many Markov models on a regular 
lattice described by the neighborhoods in Besag [4] involving coding sets or 
"unilateral" structures, there is typically no loss of generality in building 
a collection of concliques C±, . . . ,C q from such basic concliques. We illustrate 
Proposition 4.1 with some examples. 

Example 2.1 (Continued). Under the four-nearest neighbor structure 
in Z 2 , we have M = {±(0, 1)', ±(1, 0)'} = ±M, m 1 =m 2 = 1, A = diag(2,2) 
and q* =4, so there are four basic concliques {CT}j =1 determined by the 
vectors 

a! = (0,0)', a 2 = (l,l)', a 3 = (l,0)', 84 = (0,1)'. 

Because a 2 — ai + A • s = (1, 1)' + 2s ^ ±A4 for any s E Z 2 , HsH^ < 1, then 
C\ = C\ U C| is a conclique, and, similarly, so is C 2 = C| U C|. Additionally, 
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Proposition 4.1 shows also that C±, C2 cannot be further merged so that 
these represent the previously illustrated minimal conclique cover. 

Example 2.2 (Continued). Under the eight-nearest neighbor structure 
in I?, we have that M = {±(0, 1)', ±(1, 0)', ±(1, 1)', ±(1, -1)'} and the basic 
concliques {Cj}1j =1 are the same as in Example 2.1 and correspond to Be- 
sag's [4] coding sets. However, these basic concliques cannot be merged into 
larger concliques by Proposition 4.1 and hence match the minimal cover of 
four concliques as illustrated previously (i.e., Cj =Q). 

Example 4.1. Under a "simple unilateral" neighbor M = {(0, — 1)', 
(-1,0)'} in 1? (cf. [4], Section 6.2), the basic concliques are again the same 
and Proposition 4.1 gives C\ = C* U C2 , C2 = C3 U C\ as a minimal conclique 
cover. 

4.2. Asymptotic framework. We now consider a sequence of sampling 
regions R n indexed by n. For studying the large sample properties of the 
proposed GOF statistics, we adopt an "increasing domain spatial asymp- 
totic" structure [10], where the sampling region R n becomes unbounded 
as n — > 00. Let Ro be an open connected subset of (—1/2, l/2] d containing 
the origin. We regard Rq as a "prototype" of the sampling region R n . Let 
{A n }n>i be a sequence of positive numbers such that X n — > 00 as n — > 00. We 
assume that the sampling region R n = X n Ro is obtained by "inflating" the 
set Ro by the scaling factor A n (cf. [40]). Since the origin is assumed to lie 
in Rq, the shape of R n remains the same for different values of n. To avoid 
pathological cases, we assume that for any sequence of real numbers {a n } n >i 
with a n — > 0+ as n — > 00, the number of cubes of the lattice a n Z, d that inter- 
sect both Rq and Rq is 0((a n )~^ d ~^) as n — > 00. This implies that, as the 
sampling region grows, the number of observations near the boundary of R n 
is of smaller order 0(N n d 1 ^ d ) than the total number N n of observations 
in R n so that the volume of R n , N n and the number of interior locations 
are equivalent as n — > 00. The boundary condition on Rq holds for most re- 
gions Rn of practical interest, including common convex subsets of W* 1 , such 
as rectangles and ellipsoids, as well as for many nonconvex star-shaped sets 
in M. d . (Recall that a set A C M. d is called star-shaped if for any x £ A, the 
line segment joining x to the origin lies in A.) The latter class of sets may 
have a fairly irregular shape. See, for example, [38, 43] for more details. 

We want to assess the GOF of the process model specification (1.2), under 
either the simple or composite hypothesis sets of Section 3. As described in 
Section 3.1, we suppose that the spatial process is observed at locations 
on the integer grid Z d that fall in the sampling region R n producing a set 
of sampling locations Sn u (indexed by n). To simplify notation, we will 
use S n rather than the more cumbersome Sn u and S n nt rather than S 1 ^. 
Similarly, we will use Wj n to denote the empirical distribution of generalized 
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spatial residuals for the jth conclique under a simple hypothesis as given 
by (3.3) with N = N n and T± n , . . . ,T4 n , the corresponding test statistics of 
(3.4)-(3.7). Also, W jn , and ?i n ,...,T4 n will denote the quantities in (3.9) 
and (3.10) with N = N n . 

4.3. Results for the simple testing problem. For studying the asymp- 
totic distribution of the test statistics T± n , . . . ,T^ n under the null hypothe- 
sis Hq(S), we shall make use of the following condition, which imposes the 
structure on the concliques described in Section 4.1. 

Condition (C.l): Each conclique C\,...,C q is union of basic concliques 
C|,..., C*» as in (4.1). Namely, for each j = l,...,q, there exists Jj C 
{l,...,q* = det(A)} where Cj = \Ji e j j C* and the index sets {J7j}j =1 are 
disjoint. 

The following result gives the asymptotic null distribution of conclique- 
wise empirical processes W n = (W\ n , . . . , W qn )' based on the scaled and cen- 
tered empirical distributions Wj n (u), u G [0,1], as in (3.3). Note that, while 
each individual empirical process Wj n can be expected to weakly converge 
to a Brownian bridge under Hq(S) (cf. [3]), the limit law of W n will not 
similarly be distribution-free due to the dependence in observations across 
concliques. In particular, the null model F$ Q influences the asymptotic co- 
variance structure of W n . 

Let C%o denote the collection of bounded vector-valued functions f = 
(/i,..., f q )' : [0, 1] — >■ defined on the unit interval. Also, let |B| denote 
the size of a finite set BcR. 

Theorem 4.2. Suppose that condition (C.l) holds. Then, there exists 
a zero-mean vector-Gaussian process W(it) = (W\{u), . . . , W q {u))' ,u € [0, 1], 
with continuous sample paths on [0,1] (with probability 1) such that 

W n Aw as n — > oo 

as elements of Further, P(W(n) = 0) = 1 for u = 0, 1 and the q x q 
covariance matrix function of W is given by 

( det(A) 

{mm{u,v} -uv), if j = k, 

EW 3 {u)W k {v) = \ J m 

for < u, v < 1, 1 < j,k < q and 

o-ij(u,v)= Yl {P[U (0) < u, U {&i - a t + As) <v]-uv} 
sez d ,\\s\\ oc <i 

x K(ai - &i + As G ±M). 
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The indicator function !(•) above pinpoints terms in the covariance ex- 
pression which automatically vanish by the independence of residual vari- 
ables U(s) within conclique structures (Theorem 2.1). For example, when 
it is possible to combine two basic concliques C* and C£ , i 7^ /, into a larger 
conclique, Proposition 4.1 gives that, for all ||s||oo < 1, it holds that aj — a; + 
As ^ M. and so above I(aj — a^ + As € ±.M) = 0. All sums in the limiting 
covariance structure then involve only a finite number of terms. 

As a direct implication of Theorem 4.2, we get the following result on the 
asymptotic null distribution of the test statistics Ti n , . . . ,T^ n . 

Corollary 4.3. Under the conditions of Theorem 
Tj n -A ipj (W) as n — > oo 



for j = 1, . . . , 4, where the functionals ' ipj 's are defined by 
^i( f ) = ma* SU P \fj( u )l 

l<J<9 ue [0,l] 




for f = (/i, ...,f q )'e Clc, and for a given r G [l,oo). 



4.4. Results for the composite testing problem. As for the simple test- 
ing problem, here we first derive the asymptotic null distribution of the 
conclique-wise empirical processes W n = (Wi n , . . . ,W qn )' based on scaled 
and centered empirical distributions Wj n (u), u € [0,1], in (3.9). 

Note that the estimator 9 n appears in each summand in Wj n through 
the estimated generalized spatial residuals (3.8). In such situations, a com- 
mon standard approach to deriving asymptotic distributions of empirical 
processes is based on the concept of uniform asymptotic linearity in some 
local neighborhood of the true parameter value #o (cf. [36, 46]). However, 
this approach is not directly applicable here due to the form of the condi- 
tional distribution functions in (3.1) when considered as functions of 6 £ 0. 
To establish the limit distribution, we embed the empirical process of the 
estimated generalized residuals in an enlarged space, namely, the space of 
locally bounded q-dimensional vector functions on [0, 1], equipped with the 
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metric of uniform convergence on compacts, and then use a version of the 
continuous mapping theorem; the argument details are provided in [32]. 

We require some notation and conditions in addition to those introduced 
in the earlier section. Letting again \B\ denote the size of a finite set B, 
define the strong mixing coefficient of the process {F(s):seZ d } by 

a(a; b) = sup{|P(A nfl)- P(A)P{B)\: A G V(Si),B G V(S 2 ), 

\Si\<b,\S 2 \<b,d(S 1 ,S 2 )>a,S 1 ,S 2 cZ d }, 

where T>(S) = cr(Y(s) : s G S) generically denotes the c-algebra generated by 
variables Y(s) with locations in S C Z, d , d(Si,S 2 ) = inf{||s — t||i : s G Si, 
t G S 2 }, |[x||i = Yli=i \ x i\ f° r x = { x ii ■ ■ ■ i x d)' € M. d , and P(-) represents 
probabilities for the process. Write Fg (■{■) and Fg 1 ^ (-|-) to denote p X 1 
vectors of first order partial derivatives of Fg(-\-) and F ~(-\-) with respect 
to 9, when these exist. Let U [Q) = (1 - A{Q)) ■ Fg(Y(0)\{Y(t) : t G M}) + 

A{Q) ■ Fq (Y (0)\{Y {t) : t G M}), and denote uj 1} (0) € K p as the vector of 
partial derivatives of Ug(Q) with respect to 9, when this exists. 
Condition (C.2): 

(i) There exist constants 5o G (0, 1), Co G (0,oo) such that 

\P(Ue(0) < u) - P(UeM <v)\< c [\\6 - 9 \\ + \u- v\] 
for all < u,v < 1 and 9 G G satisfying max{||0 — 6q\\, \u — v\} < 5q. 

(ii) sup{||F, (1) (y|x)|| + ||F e (1) "(y|x)|| : \\9 - 9 \\ < 5 ,y G R,x G MP} < c . 

(iii) E{sup lle _ eoll<s \\U e 1 \o)-U e 1 J(0)\\} = o(5) asJ^O. 

Condition (C.3): Suppose that the joint distribution of (Uq (0),XJq (0)) is 
absolutely continuous with respect toLx/j with Radon-Nikodym derivative 
f(u,x), where L is the Lebesgue measure on R, and \i is a cr-finite measure 
on W . Suppose that 

lim sup / ||x||/(it, x) dfi{x) = 
^^ueCo.i) J[|x|[>* 

and 

y ||x|| • sup{|/(ii,x) — /(u,x)| : I it — v| < 5}d/x(x) — > 

as 5 -> 0+. 

Condition (C.4): 

(i) There exist zero-mean random variables {V(s) : s G 7j d } such that 

N l J\9 n - Oo) = iV" 1 / 2 V ( s ) + OpCO- 

sG5„ 

(ii) For each s G the variable V(s) = (T-i(s), . . . , V p (s))' is 2?(s + M)- 
measurable. 
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(iii) There exist a G (2,oo), « > such that sup{£||V(s)|| 2+K : s G Z d } < 
oo and 

oo oo 

^/-Vj;l) K/(2+K) <oo, ^/( 2r - 1 )a(i;2r-l) 1 / a < oo 

for some integer r satisfying r > (p+ 1)/(1 — a -1 ). 

(iv) £ = lim n _ s . 0O Var(iV n ^ sg5 ^ V(s)) exists and is nonsingular. 

Conditions (C.2) and (C.3) are exclusively used for handling the effects 
of the perturbation of the empirical process of the generalized residuals due 
the estimation of 9. The first displayed condition in (C.3) is an uniform in- 
tegrability condition, while the second one is a continuity condition on the 
densities /(-,■) (in u) in a weig hted L l {n) -norm. Without loss of generality, 
we shall suppose that /(n,x) = for all u ^ (0, 1) except on a set of x- values 
with /i-measure zero. Condition (C.4) allows us to relate the limit law of 
the (unperturbed) empirical process part with the variability in estimat- 
ing 9 by 9 n . If the conditional model specification is such that the spatial 
process satisfies Dobrushin's uniqueness condition (cf. [20]), then the MRF 
is strongly mixing (actually, ^-mixing) at an exponential rate and, hence, 
mixing conditions in (C.4) trivially hold. 

Theorem 4.4. Suppose that conditions (C.1)-(C.4) and the composite 
null hypothesis Hq(C) hold. Then, there exist a zero-mean vector- Gaussian 
process W(?x) = (W\{u), . . . , W q (u))' , u G [0, 1], with continuous sample paths 
on [0,1] (with probability 1) and a random variable Z = (Z±, . . . , Z p )' ~ 
N p (0,T,), both defined on a common probability space, such that as n— >-oo, 



W„4W + 1-Z' J x/(-,x)d/i(x) 



as elements of Cto, where 1 = (1, . . . , 1)' G M q . The q x q covariance matrix 
function of W is as in Theorem 4-2 and for j = 1, . . . , q, k = 1, . . . ,p and 
«G(0,1), 

The following result is a direct consequence of Theorem 4.4 and gives the 
asymptotic distribution of the test statistics under the composite null Hq(C). 

Corollary 4.5. Under the conditions of Theorem 4-4: 

fj n -A (fj + 1 • Z' J x/(-, x) dfi(x)^j as n -> oo 

for j = 1, . . . , 4, where the functionals ipj 's are as defined in (4-2)- 
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Under the composite null Hq(C), the limiting distributions involved are 
not distribution-free (i.e., depending on the true model c.d.f. Fq q in a com- 
plex covariance structure) . Empirical processes based on PIT residuals with 
parameter estimates are known to exhibit this behavior in other inference 
scenarios with time series and independent data (cf. [17]), and often two 
general approaches are considered for implementing GOF tests [37]: resam- 
pling or Khmaladze's [33] martingale transformation. The latter involves 
a type of continuous de-trending to minimize effects of parameter estimation 
and has been applied to obtain asymptotically distribution-free tests with 
other model checks using residual empirical processes based on estimated 
parameters (cf. [34, 35]). In particular, Bai [3] justified this transformation 
for tests in parametric, conditionally specified (continuous) distributions for 
time series, but considered only one empirical process of residuals. If modi- 
fied to the spatial setting, this result would entail a transformation of WjN 
from one conclique j = 1, . . . , 1 so that its limiting distribution is Brownian 
motion and distribution-free under Hq(C). The complication here is that 
with residual empirical processes from multiple concliques, after applying 
a conclique-wise transformation, the resulting limit distribution of a test 
statistic under Hq(C) would not be distribution- free due to dependence 
across concliques (akin to Theorem 4.2 in the case of no parameter esti- 
mation). Another option might be to use plug- in estimates of the covariance 
structure, using, for example, that asymptotic variances of maximum likeli- 
hood and pseudolikelihood estimators (i.e., S in Theorem 4.4) are known for 
some Markov field models [21]. But one would also have to estimate other 
complicated covariances in the limiting distribution of Theorem 4.4, which 
might be possible with subsampling variance estimation [43]. 

Spatial resampling methodologies, such as the block bootstrap (cf. [39], 
Chapter 12), might also be used to approximate sampling distributions of 
GOF statistics based on spatial residuals and knowledge of the limit distri- 
butions in Theorem 4.4 could be applied to toward justifying such bootstrap 
estimators. Simulations in Section 5 also suggest that the finite sample ver- 
sions of the GOF statistics appear to converge fairly quickly to their limits, 
at least in the case of simple null hypotheses. This implies that, in ap- 
plication, large-sample bootstrap approximations of finite-sample sampling 
distributions may be reasonable. The theoretical development of a spatial 
bootstrap for our GOF statistics is outside of the scope of this paper, but 
in Section 6 we use a parametric spatial bootstrap to calibrate GOF test 
statistics for a composite null hypothesis. 

5. Numerical results. Here we provide a small numerical verification of 
the large sample distributional results in the simple null hypothesis case, 
considering observations generated from a conditional Gaussian MRF on Z 2 
with a four-nearest neighbor structure specified by Ai = {±(0, 1)', ±(1, 0)'} 
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as in Example 2.1. The conditional model family (3.1) of Y(s) given {^(t) : 
t G J\f(s)}, s € 1? (A/"(s) = s + M), is normal with mean //^(s) = a + 
r/^ te ^ s - ) [y(t) — a], and variance r 2 > 0, where E(Y(s)) = a£R is the 
marginal process mean and \rj\ < 0.25 denotes a dependence parameter. In 
total, the model parameters 9 are (a,T,rj)'. 

5.1. Limit distributions under a simple null hypothesis. We first examine 
asymptotic null distributions of GOF test statistics in the simple testing 
problem Hq(S) : (a,r, r/)' = (ao,To, of Section 3.2 with residuals (3.2) 
given by U(s) = <3?[{Y(s) — /U ao ,»yo( s )}/' r o]- Here <&(•) denotes the standard 
normal cumulative distribution function, and, for simplicity, we will write 
hypothesized parameters chq,tq,i]q as a,r, rj in the following. 

As described in Section 3.1, the four-nearest-neighbor structure produces 
a minimal cover of two concliques C\,C2 (cf. Example 2.1), each of which 
is a union of two basic concliques C*, . . . ,C| provided in Section 4.1. These 
concliques yield an empirical distribution process W n = (Wi n ,W2n)' and 
GOF test statistics T ln , . . . ,T 4n as in (3.4)-(3.7). By Theorem 4.2, W n has 
a mean-zero Gaussian limit W = (Wi, W^)' with covariances 

(5.1) W,-(«)W,H - | g[p(Xi < $ -i (u) X2 < _ uu]j if ^ fej 

u,d£ [0, 1], j, A; G {1, 2}, where vectors (Xi, X2) in (5.1) are bivariate normal, 
with marginally standard normal distributions and correlation —rj. Hence, 
under the simple null hypothesis, the limit process depends on (a,T, r))' only 
through the dependence parameter 77, which we denote by writing W = W„. 

To understand the distribution of ^(W.^), j = 1,2,3,4, as the asymp- 
totic limit of GOF statistics Tj n under Corollary 4.3, we simulated from 
the theoretical Gaussian process as follows. For each value of rj = 
0,0.1,0.24, we generated 50,000 sequences of mean-zero bivariate Gaus- 
sian variables (Wi(i/3001), W2(z/3001)), i = 0, . . . , 3001, with covariance 
structure (5.1) over a grid in [0,1]; the sequence length of 3002 was dic- 
tated by computational stability. These provide approximate observations 
of Wjj, with rj values chosen to reflect no, weak and strong forms of positive 
spatial dependence. Cumulative distribution functions of each functional 
(pi(W„), . . . , (p^(yV„) were then approximated from W^-realizations. The 
resulting distribution curves appear in Figure 1 for rj = 0.1 and rj = 0.24, 
with </?3(Wn) and ^(W,,) computed using r = 2 in (4.2). 

5.2. Comparisons to finite sample distributions. To compare the agree- 
ment of finite sample distributions of 2Vzv under the simple null hypothesis 
with their limit distributions yjj(W^), j = 1, ... ,4, we simulated samples on 
two grid sizes, a 10 x 10 grid having N = 100 locations and a 30 x 30 grid hav- 
ing N = 900. Here, we simulated 50,000 realizations of conditional Gaussian 



18 



M. S. KAISER, S. N. LAHIRI AND D. J. NORDMAN 

Functional # 1 Functional # 2 




~i 1 1 1 1 n 1 1 1 1 r~ 

0.0 0.5 1.0 1.5 2.0 2.5 0.0 0.5 1.0 1.5 2.0 2.5 



Functional # 3 Functional # 4 




h i i i i ^ i i i i r~ 

0.0 0.5 1.0 1.5 2.0 2.5 0.0 0.5 1.0 1.5 2.0 2.5 



Fig. 1. Cumulative distribution functions i 7 ' ¥ , j (w^)(™) = ^fe(Wij) < w], w £ R, for 
limit functionate (pi (W^), . . . , ^(W,,) for r\ — 0.1 (dashed) and n — 0.24 (solid). 

samples (setting a = and r = 1 with no loss of generality) and evaluated 
functionals Tin ?4JV to approximate the finite-sample distributions of 
these GOF statistics. Figure 2 shows the difference between the quantiles of 
the limit (^(W^) and those of T2N for rj = 0.1 and r/ = 0.24; the agreement 
among quantiles for functional 2 is quite good even though this plot was 
one exhibiting the largest quantile-mismatches among the four GOF func- 
tionals. Table 1 shows the proportion of GOF statistics Tjn falling above 
the 95th and 99th quantiles of the corresponding limit </fj(Wn) distribution, 
j = 1,2,3,4. The agreement between the finite-sample and theoretical limit 
distributions is again close in Table 1. 

For various sample sizes and dependence parameters, Table 2 compares 
the finite-sample distributions of the four GOF statistics {Tjjv}|=i against 
their limiting distributions <^,-(W) in terms of a Kolmogorov-Smirnov -Dks 
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Fig. 2. Difference in quantiles for </?2(W 7) ) and T^jv when N = 100 (dashed line) and 
900 (solid line) for rj = 0.1, 0.24. Pointwise 95% confidence bands (dotted) indicate the 
Monte Carlo error in each difference. 



and a Cramer-von Mises-like -Dcm distance metric, defined by 
D KS (X,Z)=sn V \F x (t)-F z (t)\, 

1/2 

D CM (X,Z) = 



J \F x (t)-F z (t)\ 2 dt 



relative to the cumulative distributions Fx,Fz of arbitrary random vari- 
ables X, Z. To interpret the relative values of these metrics in assessing the 
distributional distance between TjN and pj(W n ), it is helpful to reference 

Table 1 

Proportion of GOF statistics TjN from a conditional Gaussian model falling above the 
95th and 99th quantiles (denoted §95,,, and q99, v ) of the their limit (pj(W v ) distribution, 
j — 1, 2, 3, 4, for sample sizes N = 100 and N — 900 and with dependence parameters 

77 = 0,0.1,0.24 



V 


N 




% of T jN 


> <795,77 






% of T jN 


> <?99,T| 




J = l 


2 


3 


4 


J = l 


2 


3 


4 





100 


4.67 


4.38 


5.09 


4.97 


0.90 


0.90 


0.98 


0.91 





900 


5.11 


4.91 


4.95 


4.86 


1.03 


1.09 


1.03 


1.03 


0.1 


100 


4.60 


4.60 


4.88 


4.92 


0.94 


0.95 


0.95 


1.08 


0.1 


900 


5.11 


5.13 


5.05 


5.08 


1.08 


1.13 


1.07 


1.15 


0.24 


100 


4.52 


4.57 


5.02 


5.06 


0.80 


0.76 


0.86 


0.92 


0.24 


900 


4.92 


4.97 


4.86 


5.03 


0.97 


0.96 


0.93 


0.95 
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Table 2 

Computed values (x 1000) from distance metrics comparing finite- sample distributions of 
statistics Tin, ■ ■ ■ , Tijv to their limiting distributions ip±(W v ), . . . , (^(W,,) 



Dks (tp, (W„ ),T jN ) Dcm (<Pi (W„ ),T jN ) 





N 


i — i 


2 


3 


4 


J = 1 


2 


3 


4 





100 


19.6 


23.0 


8.9 


9.1 


12.9 


14.8 


3.7 


3.2 





900 


6.7 


9.7 


3.4 


4.8 


3.8 


4.7 


1.0 


1.4 


0.1 


100 


24.0 


27.7 


5.3 


5.4 


16.5 


18.2 


2.2 


1.7 


0.1 


900 


9.9 


10.0 


6.5 


5.7 


5.3 


4.7 


2.5 


1.8 


0.24 


100 


21.6 


25.2 


7.2 


7.0 


14.7 


15.6 


2.7 


2.2 


0.24 


900 


9.1 


8.2 


4.2 


3.8 


4.8 


4.8 


1.4 


1.4 








),¥>j(W„ 2 


)) 


Dcm 


(<Po ( W m 


), ip j(W,. 


.)) 


Vi 




j = l 


2 


3 


4 


J = 1 


2 


3 


4 





0.1 


14.4 


11.9 


14.9 


13.4 


7.0 


5.7 


7.7 


6.1 


0.1 


0.24 


70.0 


59.6 


90.3 


72.9 


49.9 


35.0 


50.4 


33.1 





0.24 


81.8 


69.1 


102.3 


84.3 


56.6 


40.5 


58.0 


38.9 



Dks j Dqm values for comparing the distributions of tpj ( W J?1 ) and ip<j (W m ) 
over parameters r]i 7^ 7/2, which Table 2 also provides. Generally, the conver- 
gence of the finite-sample distributions Tjn to their limits ipj(W v ) appears 
to occur fairly uniformly over different dependence parameters rj and, rel- 
ative to the distributional differences among different limits [e.g., y?j(W r?1 ) 
and yj(W^ 2 )], the agreement in distributions of Tjn and ^(W^) is quite 
close even for samples of size 100. 

5.3. Power of GOF statistics under simple null hypothesis. Under the 
simple null Hq(S) : (a,T,rj)' = (0, 1,0)', we next consider the power of GOF 
tests based on statistics Tin, ■ ■ ■ ,T^n computed from conditional Gaussian 
data generated with n = 0.1 and n = 0.24 and a = 0, r = 1. This gives an 
idea of the power in testing a hypothesis of no spatial dependence, when 
the data exhibit forms of positive dependence, both fairly weak (rj = 0.1) 
and strong (n = 0.24). For a given GOF statistic Tjn from a sample of 
size N = 100 or N = 900, a size 7 test is conducted by rejecting Hq if Tjn 
exceeds the 1 — 7 quantile of the limit distribution </?(W^ = o) under the null 
hypothesis. Figure 3 plots power versus size 7 for these tests when n = 0.1 
and n = 0.24, based on 50,000 simulated data sets. Power is low under the 
alternative n = 0.1, as might be expected, but considerably higher when 
rj = 0.24. Tests with functionals T2n,T4n (based conclique-wise averages of 
GOF statistics) tend to perform similarly and exhibit slightly more power 
than tests with Tijv,T3at (based conclique-wise maxima of GOF statistics). 
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0.1 



o 
a 




o 

CD 



77 = 0.24 



jit j)M***'* M " ,, "" l,,M ' 



i 1 ' ' 



If H 



MitUiMiii'M 



0.10 

size 



Fig. 3. Plots of power versus size 7 for GOF tests of Ho : 77 = in conditional Gaussian 
models (fixed a = 0, r = 1) based on functionals Tin , ■ ■ • , TW , determined by data generated 
under 77 = 0.1, 0.24. In these power versus size curves, each functional is numbered 1~4 
under sample sizes N — 100 (grey) and N = 900 (black). 



6. An application to agricultural field trials. 

6.1. The problem. Besag and Higdon [5] present an analysis of six agri- 
cultural field trials of corn varieties conducted in North Carolina using a hi- 
erarchical model that included an intrinsic Gaussian MRF as an improper 
prior for spatial structure. An intrinsic Gaussian MRF results from fixing 
dependence parameters at the boundary of the parameter space. In discus- 
sion of this paper, Smith [44] raised the question of what diagnostics were 
available to examine potential evidence for spatial structure based on the 
available data, and presented variograms of three of the trials. Kaiser and 
Caragea [28] used data from these same three trials to illustrate a model- 
based diagnostic they called the 5-value. Questions about the spatial struc- 
ture suggested by the data included the possibilities of nonstationarity and 
directional dependencies. Here, we use data from all six trials to examine the 
question of whether a simple model with constant mean and unidirectional 
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dependence can be rejected as a plausible representation of spatial structure. 
Our question is simply one of whether a basic Gaussian MRF with constant 
mean and a single dependence parameter could be rejected as a possible 
data generating mechanism for the data, not whether it might be be most 
preferred model available. 

Each field trial consisted of observations of yield from 64 corn varieties 
with each variety replicated 3 times in each trial. The spatial layout of each 
trial was essentially that of a 11 x 18 regular lattice, although the last column 
of that lattice contained only 5 locations. After subtracting variety by trial 
means in the same manner as [28, 44], we deleted the last column to obtain 
a rectangular 11 x 17 lattice containing 187 observations for each trial. We 
assumed a four-nearest-neighborhood structure but without use of a border 
strip, so that locations had a variable number of neighboring observations, 
4 for each of the 135 interior locations, 3 for each of the 48 edge locations, 
and 2 for each of the four corner locations. 

6.2. The model. Although each trial should nominally have marginal 
mean zero, to examine a full composite setting we fit a model with con- 
ditional Gaussian distributions having expected values {/^(sj) :i = l,...,n} 
and constant conditional variance r 2 where, with N{ denoting the neighbor- 
hood of location Sfj i = 1, . . . , n, 

(6.1) fi(si) = a + i] {v( s i) ~ «}■ 

The joint distribution of this model is then Gaussian with marginal means a 
an n-vector with each element equal to a and covariance matrix (/ — C) _1 i\/ 
where I is the n x n identity matrix, M is an n x n diagonal matrix with 
all nonzero entries equal to r 2 and C = rjH with H an n x n matrix having 
element equal to 1 if locations s, and Sj are neighbors and other- 

wise. With this structure, the parameter space of r] can be determined to 
be (—0.2563,0.2563) based on eigenvalues of H (cf. [10]); this differs slightly 
from the parameter space for a lattice with four-nearest-neighborhood struc- 
ture wrapped on a torus due to the size of the lattice and the use of varying 
numbers of neighbors for edge locations. 

6.3. The GOF procedure. The model of expression (6.1) was fit to (cen- 
tered) data from each of the six trials using maximum likelihood estimation. 
Generalized spatial residuals were computed for each of the two concliques, 
one having 93 and the other 94 locations. Using the fitted models, a para- 
metric bootstrap procedure was used to arrive at p- values for each of the four 
test statistics introduced as Tj^;j = 1, ... ,4, in Section 3.3. For each fitted 
model (i.e., trial) 5000 bootstrap data sets were simulated using a Gibbs al- 
gorithm with a burn-in of 500 and spacing of 10, which appeared adequate 
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Table 3 

Estimates for conditional Gaussian models fit to data from six agricultural field trials; 
the point estimates for a for all trials differ from zero by at most 10~ 15 



Trial 


Point 




Interval 




T 2 


V 


a 


T 2 


?7 


1 


95.56 


0.2526 


(-10.21,10.40) 


(79.43,119.54) 


(0.2107,0.2544) 


2 


156.90 


0.1855 


(-3.19,3.42) 


(125.96,190.08) 


(0.0922,0.2257) 


3 


128.94 


0.2476 


(-7.66,7.76) 


(105.63,159.54) 


(0.1976,0.2533) 


4 


129.92 


0.2095 


(-3.57,3.74) 


(104.54,159.76) 


(0.1264,0.2380) 


5 


69.33 


0.2522 


(-8.29,8.23) 


(57.20,86.44) 


(0.2091,0.2543) 


() 


210.75 


0.2542 


(-20.57, 19.68) 


(175.39,268.45) 


(0.2136,0.2549) 



to result in convergence of the chain based on scale reduction factors [18] 
and eliminate dependence between successive data sets based on autocor- 
relations. Model (6.1) was fit, generalized spatial residuals produced and 
the four test statistics computed for each bootstrap data set, from which p- 
values were taken as the proportion of simulated test statistic values greater 
than those from the actual data sets. Bootstrap data sets were also used to 
produce percentile bootstrap intervals for parameters (cf. [14]). Percentile 
intervals were chosen because basic bootstrap intervals extended beyond the 
parameter space for r\ for each of the six trials. 

6.4. Results. Results of estimation are presented in Table 3. Intervals 
were computed at the 95% level and values for r\ are reported to four decimal 
places because estimates tended to be close to the upper boundary of the 
parameter space (0.2563). Overall, estimation was fairly similar for these 
six trials, which were conducted in different counties of North Carolina, 
including an indication of high variability in estimating these parameters, 
particularly a and r 2 . Estimates of r\ indicate moderate to strong spatial 
structure in each of the six trials, and estimates of r 2 indicate substantial 
local variability despite this structure. 

GOF p-values resulting from the parametric bootstrap procedure of Sec- 
tion 6.3 are presented in Table 4 for each of the four test statistics of Sec- 
tion 3.3. Overall these values provide no indication that we are able to 
dismiss model (6.1) as a plausible representation of the spatial structure 
present in these data. 

7. Conclusions. In this article we have introduced a practical method to 
assess the aptness of Markov random field models for representing spatial 
processes. This method is based on special sets of locations we have called 
concliques that partition the total set of observed locations such that gen- 
eralized spatial residuals within each conclique approximate realizations of 
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Table 4 

Parametric bootstrap p-values for the six agricultural field trials 



Trial 


2i 


T 2 


T 3 


T 4 


1 


0.8348 


0.7976 


0.7086 


0.7530 


2 


0.3844 


0.4182 


0.2132 


0.3262 


3 


0.0852 


0.1168 


0.1506 


0.1478 


4 


0.1656 


0.1084 


0.1426 


0.0972 


5 


0.2162 


0.1828 


0.1754 


0.2024 


6 


0.3502 


0.2382 


0.4642 


0.2984 



independent random variables on the unit interval. These generalized spatial 
residuals can be combined across nonindependent concliques in natural ways 
to produce GOF statistics that correspond to Gaussian empirical processes 
that have identifiable limit distributions. While those limit distributions 
can involve complex covariance structures, we have demonstrated that finite 
sample versions of the GOF statistics appear to converge rather quickly to 
their limits, at least in the case of a simple null hypothesis. This implies that, 
in an application, approximation of their limit distributions under a suitable 
null hypothesis will provide a useful reference distribution against which to 
compare the value of an observed GOF statistic. The composite hypothe- 
sis setting introduces a considerably more complicated situation than does 
the simple hypothesis setting, because limit laws involve covariances that 
cannot be easily determined either explicitly or numerically. In an appli- 
cation, resampling methods would seem to hold the greatest promise for 
approximating distributions of GOF statistics based on generalized spatial 
residuals. While developing spatial subsampling or block bootstraps (cf. [39], 
Chapter 12) for this purpose requires further investigation, the use of such 
resampling was illustrated in this article in the application to agricultural 
field trials. 

We wish to comment on a number of issues that involve the distinction 
between application of the GOF methodology developed and the production 
of theoretical results for that methodology. First is the issue of stationarity. 
There is nothing in the definition or construction of generalized spatial resid- 
uals, or GOF statistics constructed from them, that requires a stationary 
model. All that is needed is identification of a full conditional distribution 
for each location (1.1) that may then be used in (2.1), and assurance that 
a joint distribution having these conditionals exists. Assumptions of sta- 
tionarity made in this article facilitate the production of theoretical results 
needed to justify use of the methodology. Another issue is application to 
discrete cases. While the data examples given have considered continuous 
conditional models, we have applied random generalized spatial residuals to 
models formed from Winsorized Poisson conditional distributions [29] with 
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promising empirical results. Similar to questions of stationarity and dis- 
crete cases, there is nothing in the constructive methodology that requires 
a regular lattice or that each location have the same number of neighbors. 
Use of a regular lattice in this article again facilitates the demonstration of 
theoretical properties, but this is not needed to implement the procedures 
suggested. The application of Section 6 involved a regular lattice, but no 
border strip or other boundary conditions were imposed to render neighbor- 
hoods of equal size. It should certainly be anticipated that there may be edge 
effects on GOF statistics as developed here, just as there are edge effects on 
properties of estimators. How severe these effects might be in various set- 
tings, and whether the use of modified boundary conditions (e.g., [6]) could 
mitigate such effects is an issue in need of additional investigation. Essen- 
tially the same thoughts can be offered relative to potential sparseness that 
might occur in an application. Locations lacking neighbors entirely could be 
considered members of any conclique one chooses, and construction of GOF 
statistics would proceed unhindered. What the effects of varying degrees of 
sparseness are remains an open question. Of course, if no locations have 
any neighbors, then the methodology presented here reduces to statistics 
constructed on the basis of the ordinary probability integral transform for 
independent random variables. 

As with all GOF tests, the procedure based on generalized spatial residu- 
als developed in this article serves as a vehicle for assessing a selected model 
for overall adequacy, not as a vehicle for selection of the most attractive 
model in the first place. This is important in consideration of fitted models 
under the composite setting, in which we can think of estimation as hav- 
ing "optimized" a given model structure for description of a set of observed 
data. There may be two or more such structures that could be, with the 
best choice of parameter values possible, viewed as plausible data generat- 
ing mechanisms for a set of observations. This does not necessarily mean, 
however, that those different structures are equally pleasing as models for 
the problem under consideration. 

Finally, we mention a connection with the assessment of /c-step ahead 
forecasts in a time series setting. Let {Xt;t € Z} denote a series of random 
variables observed at discrete, equally spaced, points in time. The proba- 
bility integral transform with distributions conditioned on the present and 
past has been used to construct fc-step ahead residuals Ut+k = Ft+k\t(Xt+k) , 
where the conditioning in F is on {X t , X t -i, . . .} (e.g., [15, 16, 19]). While 
our use of the probability integral transform is similar to what is done in 
this time series setting, the conditioning requirements are quite distinct. In 
the spatial setting, two spatial residuals are independent only if neither is 
in the conditioning set of the other (i.e., are both in the same conclique). In 
time series fc-step ahead forecasts, two values, XJ% and Uj, will be indepen- 
dent only if either X{ is in the conditioning set of Xj, or vice versa. The 
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difference stems from the use of full conditionals in spatial Markov random 
field models, rather than the sequential conditionals in the time series con- 
text which need not invoke a Markov property at all. The approach taken 
to the development of theoretical results in this article could potentially be 
used in the time series setting, but the modifications require further inves- 
tigation. 

8. Proof of generalized spatial residual properties. As Theorem 2.1 pro- 
vides the main distributional result for generalized spatial residuals (2.1) 
from concliques, which are fundamental to the GOF test statistics of Sec- 
tions 3 and 4, we establish Theorem 2.1 here. The proofs of other results 
from the manuscript are provided in supplementary material [32]. 

Let Q denote a finite subset of conclique CcZ rf with \Q\ = I > 2, and 
let Xq = UseQ^ : 1 e -^"( s )} = { s i> • • ■ j s l}, L >1, be the finite index set of 
all neighbors of sites in Q; additionally, enumerate the I elements of Q as 
Q = {si + l, . . . , si + l}, say. With respect to the enumeration of Xq and Q, 
let -Fi(-) denote the marginal c.d.f. of Y(si), and let Fj(-), 2 < j < L + I, 
denote the conditional c.d.f. of Y(sj) given Y(si), . . . ,Y(sj-i); define the 
function Fj~(-) by the left limits of Fj(-). By the randomized PIT [7], {(1 — 
A(sj)) ■F j \yla j )] + A(8j) -Fr[Y(sj)] :j = l,.. .,L + l} are i.i.d. Uniform (0, 1) 
random variables. 

For any i, k G {L + 1, . . . , L + /}, variables Y(si) and Y(sk) belong to the 
conclique Q so that all neighbors of Y(si) and Y(sk) are among {Y^Sj)}^^. 
By the Markov property (1.2), Fj[Y(sj)) = F[Y(sj \{Y(s) : s € Af(sj)})] holds 
and we may equivalently write (2.1) as U(sj) = (1 — A(sj)) ■ Fj\Y(sj)] + 
A(sj) • FJ [y(sj)] for any j £ {L + 1, . . . ,L + I}, though these relationships 
may not necessarily hold for j = 1, . . . , L. Hence, {C/(s):s € Q} are i.i.d. 
Uniform (0,1) variables for any arbitrary finite subset Q of C. 

Acknowledgments. The authors are very grateful to an Associate Editor 
and three referees for thoughtful comments and suggestions which signifi- 
cantly clarified and improved the manuscript. 

SUPPLEMENTARY MATERIAL 

Proofs of main results for spatial GOF test statistics 

(DOI: 10.1214/11-AOS948SUPP; .pdf). A supplement [32] provides proofs of 
all asymptotic distributional results from Section 4, regarding the conclique- 
based spatial GOF test statistics in simple and composite null hypothesis 
settings (Proposition 4.1, Theorem 4.2, Corollary 4.3, Theorem 4.4, Corol- 
lary 4.5). The proof in the composite hypothesis case is particularly non- 
standard; see Section 4.4. 
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